function [theta2,deScalingX,deScalingZ,hx_1,hx_2] = ...
    Var1Threshold_GMM_closedForm_macro_spec(VarU_x2,Cov10U_x2,Cov01U_x2,...
    x2,macros,econAct,thresholdLevel,allowDeSaling,switch_h0,switch_hx)
% This function computes the closed form solution for the VAR(1) model 
% with threshold effects when the factors and macro display measurement errors

% Extending the output from the regression filter with macro variables
nm          = size(macros,1);
[nx2,T]     = size(x2);
factors     = [macros;x2];
nx          = size(factors,1);
VarU        = zeros(nx,nx,T);
Cov10U      = zeros(nx,nx,T);
Cov01U      = zeros(nx,nx,T);
for t=1:T
    VarU(nm+1:nm+nx2,nm+1:nm+nx2,t)   = VarU_x2(:,:,t);
    Cov10U(nm+1:nm+nx2,nm+1:nm+nx2,t) = Cov10U_x2(:,:,t);
    Cov01U(nm+1:nm+nx2,nm+1:nm+nx2,t) = Cov01U_x2(:,:,t);
end
% Defining regimes
regime_1    = econAct >= thresholdLevel; %a boom
regime_2    = econAct < thresholdLevel;  %a recression

% The process for the pricing factors
Xmat          = factors(:,2:T);
% The most general case is used as the benchmark
a = [regime_1(1:T-1,1)';
    regime_2(1:T-1,1)';
    factors(:,1:T-1).*repmat(regime_1(1:T-1,1)',nx,1);
    factors(:,1:T-1).*repmat(regime_2(1:T-1,1)',nx,1)];

h0RegimeOn = 1;
hxRegimeOn = 1;
VarUa         = zeros((1+hxRegimeOn)*nx+1+h0RegimeOn,(1+hxRegimeOn)*nx+1+h0RegimeOn);
Amat          = zeros(nx,(1+hxRegimeOn)*nx+1+h0RegimeOn);
for t=1:T-1
    VarUa(2+1:2+nx,2+1:2+nx)           = VarUa(2+1:2+nx,2+1:2+nx) + VarU(:,:,t)*regime_1(t,1);
    VarUa(2+nx+1:2+2*nx,2+nx+1:2+2*nx) = VarUa(2+nx+1:2+2*nx,2+nx+1:2+2*nx) + VarU(:,:,t)*regime_2(t,1);
    Amat(:,2+1:2+nx)                   = Amat(:,2+1:2+nx) + Cov10U(:,:,t+1)*regime_1(t,1);
    Amat(:,2+nx+1:2+2*nx)              = Amat(:,2+nx+1:2+2*nx) + Cov10U(:,:,t+1)*regime_2(t,1);
end

% We run the specific regressions in the VAR - equation by equation
h0_1  = zeros(nx,1);
h0_2  = zeros(nx,1);
hx_1  = zeros(nx,nx);
hx_2  = zeros(nx,nx);
Wxhat = zeros(T-1,nx);
for i=1:nx
    nn        = 1+nx+switch_h0(i)+sum(switch_hx(i,:));
    select_eq = [1 switch_h0(i) ones(1,nx) switch_hx(i,:)];
    a_eq      = zeros(nn,T-1);
    VarUa_eq  = VarUa(select_eq==1,select_eq==1);
    Amat_eq   = Amat(i,select_eq==1);
    
    idx = 1;
    % Allowing for h0_1
    a_eq(idx,:) = a(1,:);
    if switch_h0(i) == 1
        idx = idx + 1;
        a_eq(idx,:) = a(2,:);        
    else
        a_eq(idx,:) = a_eq(idx,:)+a(2,:);
    end
    % Allowing for hx_1
    a_eq(idx+1:idx+nx,:) = a(2+1:2+nx,:);
    countj = 0;
    for j=1:nx
        if switch_hx(i,j) == 1
            countj = countj + 1;
            a_eq(idx+nx+countj,:) = a(2+nx+j,:);
        else
            a_eq(idx+j,:) = a_eq(idx+j,:)+a(2+nx+j,:);
        end
    end
    
    XXinv_eq             = (a_eq*a_eq'-VarUa_eq)\eye(size(VarUa_eq,1));
    ha_eq                = (Xmat(i,:)*a_eq'-Amat_eq)*XXinv_eq;
    Wxhat(:,i)           = (Xmat(i,:)-ha_eq*a_eq)';
    
    % Adding common elements across regimes
    h0_1(i,1) = ha_eq(1,1);
    if switch_h0(i) == 1
        h0_2(i,1) = ha_eq(1,2);
    else
        h0_2(i,1) = h0_1(i,1);
    end
    hx_1(i,:) = ha_eq(1,idx+1:idx+nx);
    countj = 0;
    for j=1:nx
        if switch_hx(i,j) == 1
            countj = countj + 1;
            hx_2(i,j) = ha_eq(1,idx+nx+countj);
        else
            hx_2(i,j) = hx_1(i,j);
        end
    end
end
VarWxHat             = zeros(nx,nx);
for t=1:T-1
    omegaHat = VarU(:,:,t+1)...
        + hx_1*VarU(:,:,t)*regime_1(t,1)*hx_1' ...
        + hx_2*VarU(:,:,t)*regime_2(t,1)*hx_2' ...
        - Cov10U(:,:,t+1)*regime_1(t,1)*hx_1' - hx_1*Cov01U(:,:,t+1)*regime_1(t,1) ...
        - Cov10U(:,:,t+1)*regime_2(t,1)*hx_2' - hx_2*Cov01U(:,:,t+1)*regime_2(t,1);    
    VarWxHat = VarWxHat + 1/(T-1-2*(nx+1))*Wxhat(t,:)'*Wxhat(t,:) - 1/(T-1)*omegaHat;
end

% The process for economic activity
z                            = econAct(2:T,1)';
b                            = [ones(1,T-1); econAct(1:T-1,1)'; factors(:,1:T-1)];
varUb                        = zeros(2+nx,2+nx);
varUb(2+1:2+nx,2+1:2+nx) = sum(VarU(:,:,1:T-1),3);
ZZinv                        = (b*b'-varUb)\eye(2+nx);
hb                           = (z*b')*ZZinv;
% Compare to ordinary regression!
%[B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1) factors(:,1:T-1)'])
%STATS(1,1) gives the R2
% [B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1) factors(3:4,1:T-1)']);
% hb = [B(1,1) B(2,1) 0 0 B(3,1) B(4,1)]
%[B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1) factors(2:4,1:T-1)']);
%[B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1) factors(1,1:T-1)'.*factors(3,1:T-1)'+factors(4,1:T-1)'.*factors(4,1:T-1)' ]);
%hb = [B(1,1) B(2,1) B(3,1) 0 0 0 0]
%WzHat = R;
%[B,BINT,R,RINT,STATS] = regress(z',[ones(T-1,1) econAct(1:T-1,1)])
%STATS(1,1) gives the R2

WzHat                        = (z-hb*b)';
gama0                        = hb(1,1);
gamaz                        = hb(1,2);
gamax                        = hb(1,2+1:2+nx)';
VarWzHat                     = 0;
for t=1:T-1
    VarWzHat = VarWzHat + 1/(T-1-(nx+2))*WzHat(t,1)^2 - 1/(T-1)*gamax'*VarU(:,:,t)*gamax;
end

deScalingX = 1;
[SwxHat,checkWx] = chol(VarWxHat,'lower');
[SwzHat,checkWz] = chol(VarWzHat,'lower');
if allowDeSaling == 1
    idx = 1;
    while checkWx ~= 0 && idx <= 100
        % Reducing deScaling
        deScalingX   = deScalingX*0.0;
        for i=1:nx
            nn        = 1+nx+switch_h0(i)+sum(switch_hx(i,:));
            select_eq = [1 switch_h0(i) ones(1,nx) switch_hx(i,:)];
            a_eq      = zeros(nn,T-1);
            VarUa_eq  = VarUa(select_eq==1,select_eq==1);
            Amat_eq   = Amat(i,select_eq==1);
            
            idx = 1;
            % Allowing for h0_1
            a_eq(idx,:) = a(1,:);
            if switch_h0(i) == 1
                idx = idx + 1;
                a_eq(idx,:) = a(2,:);
            else
                a_eq(idx,:) = a_eq(idx,:)+a(2,:);
            end
            % Allowing for hx_1
            a_eq(idx+1:idx+nx,:) = a(2+1:2+nx,:);
            countj = 0;
            for j=1:nx
                if switch_hx(i,j) == 1
                    countj = countj + 1;
                    a_eq(idx+nx+countj,:) = a(2+nx+j,:);
                else
                    a_eq(idx+j,:) = a_eq(idx+j,:)+a(2+nx+j,:);
                end
            end
            
            XXinv_eq             = (a_eq*a_eq'-deScalingX*VarUa_eq)\eye(size(VarUa_eq,1));
            ha_eq                = (Xmat(i,:)*a_eq'-deScalingX*Amat_eq)*XXinv_eq;
            Wxhat(:,i)           = (Xmat(i,:)-ha_eq*a_eq)';
            
            % Adding common elements across regimes
            h0_1(i,1) = ha_eq(1,1);
            if switch_h0(i) == 1
                h0_2(i,1) = ha_eq(1,2);
            else
                h0_2(i,1) = h0_1(i,1);
            end
            hx_1(i,:) = ha_eq(1,idx+1:idx+nx);
            countj = 0;
            for j=1:nx
                if switch_hx(i,j) == 1
                    countj = countj + 1;
                    hx_2(i,j) = ha_eq(1,idx+nx+countj);
                else
                    hx_2(i,j) = hx_1(i,j);
                end
            end
        end
        
        VarWxHat    = zeros(nx,nx);
        for t=1:T-1
            omegaHat = VarU(:,:,t+1)...
                + hx_1*VarU(:,:,t)*regime_1(t,1)*hx_1' ...
                + hx_2*VarU(:,:,t)*regime_2(t,1)*hx_2' ...
                - Cov10U(:,:,t+1)*regime_1(t,1)*hx_1' - hx_1*Cov01U(:,:,t+1)*regime_1(t,1) ...
                - Cov10U(:,:,t+1)*regime_2(t,1)*hx_2' - hx_2*Cov01U(:,:,t+1)*regime_2(t,1);
            VarWxHat = VarWxHat + Wxhat(t,:)'*Wxhat(t,:)/(T-1-2*(nx+1)) - deScalingX*omegaHat/(T-1);
        end
        [SwxHat,checkWx] = chol(VarWxHat,'lower');
        idx = idx + 1;
    end
    if checkWx ~= 0
        SwxHat = NaN(nx,nx);
    end
    
    idx = 1;
    deScalingZ = 1;
    while checkWz ~= 0 && idx <= 100
         % Reducing deScaling
        deScalingZ   = deScalingZ*0.0;
        ZZinv        = (b*b'-varUb*deScalingZ)\eye(size(VarUb,1));
        hb           = (z*b')*ZZinv;
        gama0        = hb(1,1);
        gamaz        = hb(1,2);
        gamax        = hb(1,2+1:2+nx)';
        VarWzHat     = 0;
        for t=1:T-1
            VarWzHat = VarWzHat + 1/(T-1-(nx+2))*WzHat(t,1)^2 - 1/(T-1)*gamax'*deScalingZ*VarU(:,:,t)*gamax;
        end
        [SwzHat,checkWz] = chol(VarzWHat,'lower');
    end
    if checkWz ~= 0
        SwzHat = NaN;
    end
end

% Saving out
theta2 = setTheta2_threshold_macro(h0_1,h0_2,hx_1,hx_2,SwxHat,gama0,gamaz,gamax,SwzHat);

end

